Convolutions

In mathematics (and, in particular, functional analysis) convolution is a mathematical operation on two functions ($f$ and $g$) to produce a third function, that is typically viewed as a modified version of one of the original functions, giving the integral of the pointwise multiplication of the two functions as a function of the amount that one of the original functions is translated.

convolution

Convolution is similar to cross-correlation. For discrete real valued signals, they differ only in a time reversal in one of the signals. For continuous signals, the cross-correlation operator is the adjoint operator of the convolution operator.

It has applications that include probability, statistics, computer vision, natural language processing, image and signal processing, engineering, and differential equations.

The convolution can be defined for functions on groups other than Euclidean space. For example, periodic functions, such as the discrete-time Fourier transform, can be defined on a circle and convolved by periodic convolution. A discrete convolution can be defined for functions on the set of integers.

Generalizations of convolution have applications in the field of numerical analysis and numerical linear algebra, and in the design and implementation of finite impulse response filters in signal processing.

Computing the inverse of the convolution operation is known as deconvolution.

Definition

The convolution of $f$ and $g$ is written $f∗g$, using an asterisk or star. It is defined as the integral of the product of the two functions after one is reversed and shifted. As such, it is a particular kind of integral transform:

$$ \begin{align} (f * g )(t) & \, \stackrel{\mathrm{def}}{=}\ \int_{-\infty}^\infty f(\tau) g(t - \tau) \, d\tau \\ & = \int_{-\infty}^\infty f(t-\tau) g(\tau)\, d\tau. \end{align} $$

While the symbol $t$ is used above, it need not represent the time domain. But in that context, the convolution formula can be described as a weighted average of the function $f*(\tau)$ at the moment $t$ where the weighting is given by $g(−\tau)$ simply shifted by amount $t$. As $t$ changes, the weighting function emphasizes different parts of the input function.

For functions $f$, $g$ supported on only $0, \infty$ (i.e., zero for negative arguments), the integration limits can be truncated, resulting in

$$ (f * g )(t) = \int_{0}^{t} f(\tau) g(t - \tau)\, d\tau \text{ for } f, g : 0, \infty \to \mathbb{R} $$

In this case, the Laplace transform is more appropriate than the Fourier transform below and boundary terms become relevant.

Visual explanation

Example A

Visual explanation A

  1. Express each function in terms of a dummy variable $\tau.$
  2. Reflect one of the functions: $g(\tau)$→$g(-\tau).$
  3. Add a time-offset, t, which allows $g(t-\tau)$ to slide along the $\tau$-axis.
  4. Start t at −∞ and slide it all the way to +∞. Wherever the two functions intersect, find the integral of their product. In other words, compute a sliding, weighted-sum of function $f(\tau)$, where the weighting function is $g(-\tau).$

Example B

Visual explanation B

In this example, the red-colored \"pulse\", $g(\tau),$ is an even function $(\ g(-\tau)=g(\tau)\ ),$ so convolution is equivalent to correlation. A snapshot of this \"movie\" shows functions $g(t-\tau)$ and $f(\tau)$ (in blue) for some value of parameter $t,$ which is arbitrarily defined as the distance from the $\tau=0$ axis to the center of the red pulse. The amount of yellow is the area of the product $f(\tau)\cdot g(t-\tau),$ computed by the convolution/correlation integral. The movie is created by continuously changing $t$ and recomputing the integral. The result (shown in black) is a function of $t,$ but is plotted on the same axis as $\tau,$ for convenience and comparison.

Example C

Visual explanation C

In this depiction, $f(\tau)$ could represent the response of an RC circuit to a narrow pulse that occurs at $\tau=0.$ In other words, if $g(\tau)=\delta(\tau),$ the result of convolution is just $f(t).$ But when $g(\tau)$ is the wider pulse (in red), the response is a \"smeared\" version of $f(t).$ It begins at $t=-0.5,$ because we defined $t$ as the distance from the $\tau=0$ axis to the center of the wide pulse (instead of the leading edge).

In [1]:
# Bring in python image analysis libraries
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
from skimage import color
import skimage.filters as filters
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage import feature
from skimage import morphology
from skimage.draw import circle_perimeter
from skimage import img_as_float, img_as_ubyte
from skimage import segmentation as seg
from skimage.morphology import watershed
from scipy import ndimage as nd
from scipy.ndimage import convolve
from skimage import feature
import glob # for bulk file import

# Set defaults
plt.rcParams['image.cmap'] = 'gray' # Display grayscale images in... grayscale.
plt.rcParams['image.interpolation'] = 'none' # Use nearest-neighbour
plt.rcParams['figure.figsize'] = 10, 10

# Import test images
imgpaths = glob.glob("./images/*.jpg") + glob.glob("./images/*.png")
# imgpaths = glob.glob("img/*.jpg") + glob.glob("img/*.png")  Windows
# Windows has different relative paths than Mac/Unix
imgset = [mpimg.imread(x) for x in imgpaths]

# Display thumbnails of the images to make sure they loaded
plt.figure()
for i,img in enumerate(imgset):
    plt.subplot(1, len(imgset), i+1)
    plt.imshow(img, cmap = 'gray')
In [2]:
for i,img in enumerate(imgset):
    plt.figure()
    plt.imshow(img)

IOPub data rate exceeded.
The notebook server will temporarily stop sending output
to the client in order to avoid crashing it.
To change this limit, set the config variable
`--NotebookApp.iopub_data_rate_limit`.
In [3]:
for i,img in enumerate(imgset):
    img_bw = img_as_float(color.rgb2grey(img))
    plt.figure()
    plt.subplot(1, 2, 1) # One row, two columns, image 1
    plt.imshow(img)
    plt.subplot(1, 2, 2)  # One row, two columns, image 2
    plt.imshow(img_bw)

Simple Edge Detection

A simple edge detection algorithm would look for edges where there is the greatest difference amongst pixels and their neighbors. Let's try it.

In [4]:
# Find horizontal edges using a simple shifting method
def find_horizontal_edges(img):
    imgbw = img_as_float(color.rgb2grey(img))
    return np.abs(imgbw[:, 1:] - imgbw[:, :-1])

# Find vertical edges using a simple shifting method
def find_vertical_edges(img):
    imgbw = img_as_float(color.rgb2grey(img))
    return np.abs(imgbw[1:, :] - imgbw[:-1, :])
In [5]:
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 3, 1)
    plt.title('Original')
    plt.imshow(img)
    plt.subplot(1, 3, 2)
    plt.title('Horizontal Edges')
    plt.imshow(find_horizontal_edges(img))
    plt.subplot(1, 3, 3)
    plt.title('Vertical Edges')
    plt.imshow(find_vertical_edges(img))

Downsample

In [6]:
# Downsample an image by skipping indicies
def downsample_image(img, skip):
     return img[::skip,::skip]
In [7]:
# Apply to image set
for i,img in enumerate(imgset):
    img = downsample_image(img, 11) # downsample    
    plt.figure()
    plt.subplot(1, 3, 1)
    plt.title('Original')
    plt.imshow(img)
In [8]:
# Apply to image set
for i,img in enumerate(imgset):
    img = downsample_image(img, 5) # downsample    
    plt.figure()
    plt.subplot(1, 3, 1)
    plt.title('Original')
    plt.imshow(img)
    plt.subplot(1, 3, 2)
    plt.title('Horizontal Edges')
    plt.imshow(find_horizontal_edges(img))
    plt.subplot(1, 3, 3)
    plt.title('Vertical Edges')
    plt.imshow(find_vertical_edges(img))

Handmade Convolutions

As Wikipedia says (https://en.wikipedia.org/wiki/Convolution ) a convolution is a mathematical operation on two functions (f and g); it produces a third function, that is typically viewed as a modified version of one of the original functions, giving the integral of the pointwise multiplication of the two functions as a function of the amount that one of the original functions is translated.

Create some arbitrary kernels

In [9]:
def generate_kernel_a(size):
    kernel = np.zeros((size,size))
    for x in range(size):
        if x%3==0:
          kernel[x][x] = 1.0
          kernel[size-3][x] = 1.0
          kernel[size-1-x][x] = 1.0            
    kernel = kernel / float(sum(sum(kernel)))
    return kernel

kernel_a = generate_kernel_a(33)
plt.imshow(kernel_a)
Out[9]:
<matplotlib.image.AxesImage at 0x10523c358>
In [10]:
def generate_kernel_b(size):
    kernel = np.zeros((size,size))
    for x in range(size):
        kernel[x][x] = 1.0
        kernel[size-1-x][x] = 1.0
    kernel = kernel / float(sum(sum(kernel)))
    return kernel

kernel_b = generate_kernel_b(55)
plt.imshow(kernel_b)
Out[10]:
<matplotlib.image.AxesImage at 0x115674630>
In [11]:
def generate_kernel_c(size):
    kernel = np.zeros((size,size))              
    return kernel

kernel_c = generate_kernel_c(33)
plt.imshow(kernel_c)
Out[11]:
<matplotlib.image.AxesImage at 0x111fa9d30>
In [12]:
def generate_kernel_d(size):
    kernel = np.zeros((size,size))
    for x in range(size):
        for y in range(size):
          if (x+y)%2==0:
            kernel[x][y] = 1.0
    kernel = kernel / float(sum(sum(kernel)))            
    return kernel

kernel_d = generate_kernel_d(33)
plt.imshow(kernel_d)
Out[12]:
<matplotlib.image.AxesImage at 0x115951358>

Convolution - convolution is a mathematical operation on two functions (f and g)

In [13]:
for i,img in enumerate(imgset):
    imgbw = color.rgb2grey(img)
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(imgbw)
    plt.subplot(1, 2, 2)
    plt.imshow(convolve(imgbw, kernel_a))
In [14]:
for i,img in enumerate(imgset):
    imgbw = color.rgb2grey(img)
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(imgbw)
    plt.subplot(1, 2, 2)
    plt.imshow(convolve(imgbw, kernel_b))
In [15]:
for i,img in enumerate(imgset):
    imgbw = color.rgb2grey(img)
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(imgbw)
    plt.subplot(1, 2, 2)
    plt.imshow(convolve(imgbw, kernel_c))
In [16]:
for i,img in enumerate(imgset):
    imgbw = color.rgb2grey(img)
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(imgbw)
    plt.subplot(1, 2, 2)
    plt.imshow(convolve(imgbw, kernel_d))

Gaussian Blur

In [17]:
# Apply to image set
sigma = 1.0
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    plt.imshow(filters.gaussian(img, sigma,mode='reflect', multichannel=False))
In [18]:
# Apply to image set
sigma = 5.0
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    plt.imshow(filters.gaussian(img, sigma,mode='reflect', multichannel=False))
In [19]:
# Apply to image set
sigma = 10.0
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    plt.imshow(filters.gaussian(img, sigma,mode='reflect', multichannel=False))

Sobel operator

The Sobel operator (https://en.wikipedia.org/wiki/Sobel_operator) uses two 3×3 kernels which are convolved with the original image to calculate approximations of the derivatives – one for horizontal changes, and one for vertical. If we define A as the source image, and Gx and Gy are two images which at each point contain the horizontal and vertical derivative approximations respectively, the computations are as follows

Gx = np.array([[-1, 0, +1],[-2, 0, +2],[-1, 0, +1]])
Gy = -Gx.transpose()

Since the Sobel kernels can be decomposed as the products of an averaging and a differentiation kernel, they compute the gradient with smoothing.

The Sobel operator, sometimes called the Sobel--Feldman operator or Sobel filter, is used in image processing and computer vision, particularly within edge detection algorithms where it creates an image emphasising edges. It is named after Irwin Sobel and Gary Feldman, colleagues at the Stanford Artificial Intelligence Laboratory (SAIL). Technically, it is a discrete differentiation operator, computing an approximation of the gradient of the image intensity function. At each point in the image, the result of the Sobel--Feldman operator is either the corresponding gradient vector or the norm of this vector. The Sobel--Feldman operator is based on convolving the image with a small, separable, and integer-valued filter in the horizontal and vertical directions and is therefore relatively inexpensive in terms of computations. On the other hand, the gradient approximation that it produces is relatively crude, in particular for high-frequency variations in the image.

Formulation

The operator uses two 3×3 kernels which are convolved with the original image to calculate approximations of the derivatives -- one for horizontal changes, and one for vertical. If we define A as the source image, and G~x~ and G~y~ are two images which at each point contain the horizontal and vertical derivative approximations respectively, the computations are as follows:

$$\mathbf{G}_x = \begin{bmatrix} +1 & 0 & -1 \\ +2 & 0 & -2 \\ +1 & 0 & -1 \end{bmatrix} * \mathbf{A} \quad \mbox{and} \quad \mathbf{G}_y = \begin{bmatrix} +1 & +2 & +1\\ 0 & 0 & 0 \\ -1 & -2 & -1 \end{bmatrix} * \mathbf{A}$$

where $*$ here denotes the 2-dimensional signal processing convolution operation.

Since the Sobel kernels can be decomposed as the products of an averaging and a differentiation kernel, they compute the gradient with smoothing. For example, $\mathbf{G_x}$ can be written as

$$\begin{bmatrix} +1 & 0 & -1 \\ +2 & 0 & -2 \\ +1 & 0 & -1 \end{bmatrix} = \begin{bmatrix} 1\\ 2\\ 1 \end{bmatrix} \begin{bmatrix} +1 & 0 & -1 \end{bmatrix}$$

The x-coordinate is defined here as increasing in the "right"-direction, and the y-coordinate is defined as increasing in the "down"-direction. At each point in the image, the resulting gradient approximations can be combined to give the gradient magnitude, using:

$$\mathbf{G} = \sqrt{ {\mathbf{G}_x}^2 + {\mathbf{G}_y}^2 }$$

Using this information, we can also calculate the gradient\'s direction:

$$\mathbf{\Theta} = \operatorname{atan}\left({ \mathbf{G}_y \over \mathbf{G}_x }\right)$$

where, for example, Θ is 0 for a vertical edge which is lighter on the right side.

In [20]:
# Apply a Sobel filter to the image.
def roll_your_own_sobel(img):
    sobel_kernel = np.array([[-1, 0, +1],[-2, 0, +2],[-1, 0, +1]])
    Gx = convolve(img, sobel_kernel)
    Gy = convolve(img, -sobel_kernel.transpose())
    return np.sqrt(np.square(Gx) + np.square(Gy))

# This should do the same thing as skimage.filters.sobel(img)
In [ ]:
# Apply to image set
for i,img in enumerate(imgset):
    imgbw = img_as_float(color.rgb2grey(img))
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(imgbw)
    plt.subplot(1, 2, 2)
    plt.imshow(roll_your_own_sobel(imgbw))
In [ ]:
# Apply to image set
for i,img in enumerate(imgset):
    imgbw = img_as_float(color.rgb2grey(img))
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(imgbw)
    plt.subplot(1, 2, 2)
    plt.imshow(filters.sobel(imgbw))

Sepia filter

In [ ]:
#create a kernel which will be used to the original image.
sepia_filter = np.array([[.393, .769, .189],
                         [.349, .686, .168],
                         [.272, .534, .131]])

def sepia(image, filter):
    sepia_img = image.dot(filter.T) #Dot product the image and array
    #T is the transpose of the array.
    #Since our filter lines do not have unit sum, so we need to rescale
    sepia_img /= sepia_img.max()
    return sepia_img
In [ ]:
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,sepia_filter))
    except:
      pass
In [ ]:
purple_filter = np.array([[.03, .769, .189],
                         [.349, .36, .1],
                         [.55, .534, .7]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,purple_filter))
    except:
      pass
In [ ]:
filter_a = np.array([[1.0, 0.0, 0.0],
                         [1.0, 0.0, 0.0],
                         [1.0, 0.0, 0.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,purple_filter))
    except:
      pass
In [ ]:
filter_all = np.array([[1.0, 1.0, 1.0],
                         [1.0, 1.0, 1.0],
                         [1.0, 1.0, 1.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_all))
    except:
      pass
In [ ]:
filter_b = np.array([[1.0, 1.0, 1.0],
                         [0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_b))
    except:
      pass
In [ ]:
filter_c = np.array([[1.0, 0.0, 0.0],
                         [0.0, 1.0, 0.0],
                         [0.0, 0.0, 1.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_c))
    except:
      pass
In [ ]:
filter_d = np.array([[0.5, 0.5, 0.5],
                         [0.5, 0.5, 0.5],
                         [0.5, 0.5, 0.5]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_d))
    except:
      pass
In [ ]:
filter_e = np.array([[0.0, 0.0, 0.0],
                         [1.0, 1.0, 1.0],
                         [0.0, 0.0, 0.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_e))
    except:
      pass
In [ ]:
filter_f = np.array([[0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0],
                         [1.0, 1.0, 1.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_f))
    except:
      pass
In [ ]:
filter_g = np.array([[1.0, 1.0, 1.0],
                         [0.5, 0.5, 0.5],
                         [0.0, 0.0, 0.0]])
# Apply to image set
for i,img in enumerate(imgset):
    plt.figure()
    plt.subplot(1, 2, 1)
    plt.imshow(img)
    plt.subplot(1, 2, 2)
    try:
      plt.imshow(sepia(img,filter_g))
    except:
      pass

Last update October 3, 2017

The text is released under the CC-BY-NC-ND license, and code is released under the MIT license.